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ABSTRACT 

A simple numerical procedure for estimating the 
stochastic robustness of a linear, time-invariant system is 
described. Monte Carlo evaluation of the system’s 
eigenvalues allows the probability of instability and the 
related stochastic root locus to be estimated. This 
definition of robustness is an alternative to existing 
deterministic definitions that address both structured 
and unstructured parameter variations directly. This 
analysis approach treats not only Gaussian parameter 
uncertainties but non-Gaussian cases, including uncertain- 
but-bounded variations. Trivial extensions of the 
procedure admit alternate discriminants to be considered. 
Thus, the probabilities that stipulated degrees of 
instability will be exceeded or that closed-loop roots will 
leave desirable regions also can be estimated. Results 
are particularly amenable to graphical presentation. 

INTRODUCTION 

Control system robustness is defined as the ability to 
maintain satisfactory stability or performance 
characteristics in the presence of all conceivable system 
parameter variations. While assured robustness may be 
viewed as an alternative to gain adaptation or scheduling 
to accommodate known parameter variations, more often 
it is seen as protection against uncertainties in plant 
specification. Consequently, a statistical description of 
control system robustness is consistent with what may be 
known about the structure and parameters of the plant's 
dynamic model. 

Guaranteeing robustness has long been a design 
objective of control system analysis, although in most 
instances, insensitivity to parameter variations has been 
treated as a deterministic problem (see Ref. 1 for a 
comprehensive presentation of both classical and modern 
robust control). Stability (gain and phase) margins are 
useful concepts for designing robust single-input/single- 
output systems, addressing disturbance rejection and 
other performance goals in the process, and they are 
amenable to the manual graphical procedures that 
preceded the widespread use of computers. With the help 
of these computers, singular-value analysis has extended 
the frequency-domain approach to multi -input/multi- output 
systems (e.g., [2,3]); however, guaranteed-stability- 
bound estimates often are unduly conservative, and the 
relationship to parameter variations in the physical 
system is weak. Structured- singular-value analysis [4] 
reduces this conservatism somewhat, and alternate 
treatments of structured parameter variations have been 
proposed (e.g., [5-7]), though these approaches remain 
deterministic. Elements of stochastic stability [8] have 
application to robustness but have yet to be presented in 
that context. 
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The notion of probability of instability , which is 
central to the analysis of stochastic robustness, was 
introduced in Ref. 9, with application to the robustness of 
the Space Shuttle's flight control system, and it is further 
described in Ref. 10. This method determines the 
stochastic robustness of a linear, time-invariant system by 
the probability distributions of closed-loop eigenvalues, 
given the statistics of the variable parameters in the 
plant's dynamic model. The probability that all of these 
eigenvalues lie in the open left-half s plane is the scalar 
measure of robustness. 

With the advent of fast graphics workstations and 
supercomputers, the stochastic robustness of a system is 
easily computed by Monte Carlo simulation, and results 
can be displayed pictorially, providing insight into 
otherwise hidden robustness properties of the system. 
The method is computationally simple, requiring only 
matrix manipulation and eigenvalue computation, and it is 
inherently non -conservative, given a large enough sample 
space. Furthermore, the analysis of stochastic 
robustness is a logical adjunct to parameter- space control 
design methods [11-14]. Details of the approach and 
examples are given in the sequel. 

PROBABILITY OF INSTABILITY 

Consider a linear, time-invariant (LTI) system of the 
form, 

x(0 = F(p)x(0 + G(p)u(f) (1 

y(r) = H(p)x(r) (2 

where x(r), u(f), y(f), and p are state, control, output, and 
parameter vectors of dimension n y m y q , and r, 
respectively, accompanied by conformable dynamic, 
control, and output matrices F, G, and H, which may be 
arbitrary functions of p. The plant is subject to LTI 
control, 

u(r) = uc(0 - CH(p)x(r) (3 

Uc(0 is a command input vector, and, for simplicity, the (m 
x n) control gain matrix C is assumed to be known without 
error. The n eigenvalues, Xi = Oj + jco /, i = 1 to n y of the 
matrix [F(p) - G(p)CH(p)] determine closed-loop 
stability and can be determined as the roots of the 
determinant equation, 

l5l„-[F(p)-G(p)CH(p)]l = 0 (4 

where s is a complex operator and \ n is the (n x n) 
identity matrix. System stability requires that no 
eigenvalues have positive real parts. The relationship 
between parameters and eigenvalues is complicated. 
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Even if F, G, and H are linear functions of p, the 
associations between p and the X( are nonlinear, and 
there is the further possibility of products of parameters in 
the feedback term. Consequently, while small-parameter- 
variation sensitivities of the eigenvalues can be estimated 
by linear methods [10], the hopes for generally applicable 
analytic expressions are slim. 

Putting aside the mathematical intricacies, note that 
the probability of stability plus the probability of instability 
is one: 

Pr(stability) + Pr(instability) = 1 (5 

Stochastic robustness is achieved when the probability of 
stability (instability) is large (small). Since stability 
requires all the roots to be in the open-left-half 5 plane, 
while instability results from even a single right-half s 
plane root, we may write 

Pr(instability) = P = 1- jpr(o)do (6 

“ oo 

where o is an /j-vector of the real parts of the system’s 
eigenvalues, pr(a) is the joint probability density function 
of o, and the integral that defines the probability of 
stability is evaluated over the space of individual 
components of a. 

Estimating the probability of stability of a closed-loop 
system from repeated eigenvalue calculation is a 
straightforward task. Denoting the probability density 
function of p as pr(p), eq. 4 is evaluated J times with each 
element of py ,j = 1 to 7, specified by a random-number 
generator whose individual outputs are shaped by pr(p). 
This Monte Carlo evaluation of the probability of stability 
becomes increasingly precise as J becomes large. Then, 

P . . M^max^O) 

Jpr(o)do = lim J (7 

-oo J -> OO 

N{.) is the number of cases for which all elements of a are 
less than or equal to zero, that is, for which o ma x ^ 0> 
where Omax is the maximum real eigenvalue component 
in a. An important feature of this definition is that it does 
not depend on the eigenvalues and eigenvectors retaining 
fixed structures. As parameters change, complex roots 
may coalesce to become real roots (or the reverse), and 
modes may exchange relative frequencies. The only 
matter for concern is whether or not all real parts of the 
eigenvalues remain in the left-half s plane. The stable 
space of o is a hypercube with one comer at the origin and 
all other comers at various infinite points. 

There is, of course, no limitation on admissible 
specifications for the multivariate pr(p): it may be 
Gaussian or non-Gaussian, as appropriate. Rayleigh, 
correlated, and any other well-posed distributions are 
admissible, the principal challenge being to properly shape 
(and correlate) the outputs of the random-number 


generator. In practice, system parameter uncertainties 
are most likely to be bounded, as typical quality control 
procedures eliminate out-of-tolerance devices, and there 
are physical limitations on component size, weight, shape, 
etc. The rectangular (uniform) distribution is particularly 
interesting, as it readily models bounded uncertainty, and 
it is the default distribution of most algorithms for random- 
number generation. Given binary distributions for each 
parameter, in which the elements of p take maximum or 
minimum values with equal probability, the Monte Carlo 
evaluation reduces to 2 r deterministic evaluations, the 
result is exact, and the probability associated with each 
possible value of p is l/2 r . Similarly, the distribution for r 
parameters, each of which takes w values (i.e., for 
quantized rectangular distributions), can be obtained from 
w r evaluations; the probability of acquiring each value of p 
(for equally probable parameter values) is l/w r . 


Histograms and cumulative distributions for varying 
degrees of stability are readily given by the Monte Carlo 


estimate of 



, where X represents a maximum 


real eigenvalue component, and -«> < X 


histogram is a plot of 


N[(L - A) < Omax - E] 


< oo. The 
vs. X; A is 


an increment in X, N[.] is the number of cases whose 
maximum real eigenvalue components lie in the increment, 
and J is the total number of evaluations. The histogram 
estimates the stability probability density function , pr(X), 
which is obtained in the limit for a continuous distribution 
of X as A -> 0 and J -> «>. The cumulative probability 
distribution of stability , Pr(X), is similarly estimated and 


presented as 


A(Omax - £) 
J 


vs. X, the exact distribution 


being achieved in the limit as J -> oo. Consequently, 


P = 1 - Pr(0). 


(8 


The regions of varying stability degree are hypercubes in 
o space, each with one comer at the n-vector X = [I EE... 
X]T and all remaining corners at appropriate infinite 
locations. 


When has stochastic robustness been achieved? The 
answer is problem-dependent. In some applications 
involving bounded parameters, it will be possible to 
choose C such that P = 0, and that is a desirable goal; 
however, if admissible parameter variations are 
unbounded, if C is constrained, or if the rank of CH is less 
than n y the minimum P may be greater than zero. C then 
must be chosen to satisfy performance goals and one of 
two robustness criteria: minimum P, or P small enough to 
meet a reliability specification (e.g., one chance of 
instability in some large number of realizations). 

STOCHASTIC ROOT LOCUS 

While it is not necessary to plot the eigenvalues (or 
roots) of eq. 4 to determine or portray stochastic 
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robustness, stochastic root loci provide insight regarding 
the effects of parameter uncertainty on system stability. 
Consider, for example, a classical second-order system 
whose roots are solutions to the equation 

s 2 + 2£a> n s + CD n 2 = 0 (9 

Suppose that the damping ratio (£) and natural frequency 
(con) are nominally 0.707 and 1, respectively, and that 
each may be a Gaussian-distributed random variable with 
standard deviation of 0.2. Allowing first £ to vary, then 
co n , 100-sample scatter plots of the roots are obtained 
(Fig. 1). These root loci are immediately recognized as 
following the classical configurations of root locus 
construction [15], with the heaviest density of roots in the 
vicinities of the nominal values. The density of roots 
depicts the likelihood that eigenvalues vary from their 
nominal values if either damping ratio or natural frequency 
is uncertain. These stochastic root loci include branches 
on the real axis and in the right-half s plane for large 
enough variations of £ and co n . 
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Figure 1. Stochastic root loci of a second-order system with Gaussian 
damping ratio or natural frequency. £o = 0.707, o>no = 1; 100 Monte 
Carlo evaluations. 

a) Effect of C variation with 0.2 standard deviation. 

b) Effect of (On variation with 0.2 standard deviation. 

If both £ and co n are uncertain and uncorrelated (i.e., p ' 
= [£ con]T), the scatter plots become "clouds" surrounding 
the nominal values; Fig. 2a is one representation of the 
resulting stochastic root locus based on the calculation of 
4,000 samples. Further understanding can be gained by 
plotting the density of roots in a third dimension above the 
root locus plot. This is done in two steps. The first step 
is to divide the s plane into subspaces (or "bins”), as in 
Fig. 2b, and to count the number of roots in each bin as a 
sampled estimate of the root density p. The result is a 
multivariate histogram, with c and co serving as 
independent variables. Complex root bins are elemental 
areas, for which pA is defined in units of roots/unit area. 
Real root bins are confined to the real axis; hence, pL 
measures roots/unit length. 

The second step is to portray the root density 
distribution. This can be done by brightening or darkening 
the bin outlines (Fig. 2b), graphing contours of equal root 
density on the two-dimensional plot, or by plotting an 
oblique view of the three-dimensional histogram or root 
density surface, as in Fig. 2c. The plotted surfaces would 
become smoother as the number of evaluations increased. 


Numerical smoothing also can be applied (judiciously) to 
account for sampling effects on the plotted surface. For 
this paper, the graphical presentations are relatively 
crude, but it is apparent that more sophisticated graphical 
processing, including the use of false color, hidden-line 
removal, surface generation, and shading can be applied to 
good effect. Root densities along the real axis present a 
special problem for 2-D presentation, in that their 
distributions are linear, not areal; oblique 3-D views 
provide a satisfactory alternative. 



Figure 2. Stochastic root loci of a second-order system with Gaussian 
damping ratio and natural frequency. = 0.707, con 0 = 1; 4,000 Monte 
Carlo evaluations. 

a) Scatter plot. 

b) 2-dimensional binned representation. 

c) Oblique 3-dimensional representation. 

There is an ostensible relationship between pr(£) and 
p; however, the relationship may be multivalued and 
ambiguous. When considering instability, distinction must 
be made between the number of cases with right-half- 
plane roots and the number of roots in the right-half plane. 
For example, a third-order system with a complex pair of 
roots and a real root can be unstable with I, 2, or 3 roots 
in the right-half plane, yet N would be incremented by one 
in each case. A high-order system with real roots could 
be unstable with one or more roots in the same right-half- 
plane bin. Again, N would be incremented by one, 
although the bin's p depends upon the number of roots it 
contains. 

A STOCHASTIC ROBUSTNESS EXAMPLE 

Reference 16 provides a linear-quadratic-Gaussian 
(LQG) design problem with a closed-loop system that is 
nominally stable, but whose stability margins become 
vanishingly small as control and estimation gains become 
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large. That example is used here for a demonstration of 
stochastic stability robustness. An unstable second-order 
plant, 
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is to be stabilized by an LQG regulator with controller 
cost function matrices, 


Q = Q 



, R = R = 1 


(12,13 


and disturbance and measurement-error spectral density 
matrices 


W = W 



N = N = 1 


(14,15 


The corresponding LQG control and estimator gains, C 
and K, are [16] 

C = (2 + V4Tq) [1 1] = [c c] (16 

K = (2 + V4TW)[1 1] T = fk k] T (17 


If the actual control effect matrix is G = [0 |i]T rather than 
[0 1] T , eq. 4 can be expressed for this problem (with the 
state consisting of the original state and its estimate, x T 
= [xi X2 xi x2]), as, 


(5-1) -10 0 

0 (5- 1 ) (IC |iC 

-k 0 (5-1+k) -1 

-k 0 (c+k) (5-1+c) 


= 54+c353+C252+ci5+co - 0 
(18 


Using Routh’s criterion, Doyle showed that (I remaining in 
(a,b) = {(1 + 1/ck), [1 - (k + c - 4)/2ck]} is a necessary 
condition for stability to be retained [17]. 


Consider two cases with different LQG gains. In 
Case 1, c = k = 4 (the limiting case as Q and W approach 
zero), and in Case 2, c = k = 100. Because 

ci = k + c - 4 + 2(|i-l)ck, co = 1 + (l-ji)ck (19,20 

the characteristic polynomial can be expressed as 

5 4 + C35^ + C25 2 + (k+c-4-2ck)5 + (1+ck) + ^xck(25-l) = 0 

(21 

which is a root-locus problem with pck taken as the gain. 
The nominal roots are found with fi = 1, and they are 


Case 1 

*1-4 = -1,-1, -1,-1 

(22 

Case 2 

X-i-4 = -0.01,-0.01,-98,-98 

(23 


Three features are immediately evident. The root locus 
gain is proportional to ck; hence, |i has a greater effect on 
the root locus in the latter case. There is a transmission 
zero at +1/2 that will draw one root into the right half 
plane. The excess of poles over zeros is three, indicating 
that additional instability must occur for large magnitudes 
of (|i - 1). There will be either one or two unstable roots 
among those going to infinity, depending on the sign of 

(H-l). 

The stochastic root locus plots based on 3,500 Monte 
Carlo evaluations with p = | x corroborate these 
predictions (Fig. 3). It is assumed that |i is a Gaussian 
random variable with mean equalling (a + b)/2, 
representing a bias from the nominal } i used to determine 
the gains, and standard deviation of (b - a)/2. In both 
cases, the root distributions are skewed and/or multi- 
modal, and each of the branches has a pronounced peak. 
Few roots lie near breakaway points, but rather 
accumulate nearer to the transmission zero or infinity. 
Figure 3a shows three of the five possible unstable 
branches, while for the higher gain, only two branches 
reach instability. Figure 4 indicates that the resulting 
Pr(£) are non-Gaussian. The corresponding probabilities 
of instability P are 0.48 and 0.33, indicating that the 
resulting distributions are dissimilar, even though the 
standard deviations were equally scaled for each case. 
(When the Case 1 standard deviation is used with Case 
2’s gains, P climbs to 0.96.) 
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Figure 4. Histograms and cumulative probability distributions for the 
Doyle LQG counterexample. Gaussian parameter uncertainty with 
mean = (a + b)/2 and standard deviation = (a - b)/2. 
a) c = k = 4, b) c = k = 100 


Now consider two similar cases in which p is a 
random variable with uniform probability in (a,b). For 
Case 1, Figures 5 and 6 illustrate how the stochastic root 
locus and probability distributions are bounded in 
comparison with Fig. 3a and 4a. In this example, the 
bounds given by Routh's criterion are not the actual 
stability bounds, and the probability of instability is non- 
zero. For Cases 1 and 2, the probabilities of instability P 
decrease to 0.27 and 0.01, respectively, as some unstable 
values associated with the tails of the p distribution have 
been eliminated. Naturally, if p had been uniformly 
distributed just inside the actual stability boundaries 
(0.9243 < p < 1.0625, for c = 4), P would be zero. 



Figure 5. Stochastic root locus for the Doyle LQG counterexample, c = 
k = 4. Parameter uniformly distributed in (a,b) = (0.875,1.0625). 




Figure 6. Histogram and cumulative probability distribution for the 
Doyle LQG counterexample, c = k = 4. Parameter uniformly 
distributed in (a,b) = (0.875, 1.0625) 


Using Loop Transfer Recovery (LQG/LTR) [17], 
linear-quadratic (LQ) robustness can be fully recovered. 
Recovery as a function of design parameter v (W=vW 0 ) 
for Case 1 is illustrated in Fig. 7 through both singular- 
value analysis and the stochastic root locus. The LQG 
return difference function in this case is a scalar, and the 
singular value is identically the return difference function: 

l+a(s) = I+C[sI-(F-GC-KH)]'lKH[5l-F]-lG (24 

The original LQ stability margins are not fully recovered 
until v > 10,000 (Fig 7a). Figure 7b illustrates the 
mechanism of recovery: increasing v pushes two 
eigenvalues to higher frequencies and decreases the 
variation in the two roots near the origin. Based upon 
3,500 evaluations and a Gaussian p variation, the present 
analysis estimates P to be zero when v >100. 




COMPUTATIONAL ISSUES 


The validity of the Monte Carlo analysis is dependent 
on the number of eigenvalues computed, the number of 
varying parameters, their probability distributions, and 
required confidence levels. The number of evaluations 
required can be related to the number of varying 
parameters by considering uniform probability 
distributions. Quantized uniform distributions 
approximate continuous uniform distributions, approaching 
them in the limit as the numbers of discrete parameter 
values go to infinity. Given n Monte Carlo evaluations of 
a system with r continuous uniform parameters, the result 
is, at best, equivalent to results computed 
deterministically for a system with r uniform parameters 
quantized in w levels, where w -n^ r . Conversely, the 
number of evaluations should be of 0( awQ, where w is an 
acceptable level of parameter quantization , and a » 1. 
Note that in a 10-parameter case, direct equivalence to 
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ten parameter quantization levels requires over 10 billion 
evaluations, while 10,000 evaluations yield results that 
are equivalent to a quantization level less than three. 
Work remains to be done in associating small-sample 
evaluations with confidence levels of the histograms and 
the resulting probability of instability. 

If Omax is monotonic in the individual elements of p, 
then evaluation results for the binary probability 
distribution denoted by (PmimPmax) circumscribe results 
obtained for continuous or quantized distributions with the 
same limits. In this case, a conservative estimate of P is 
provided by the associated 2 r deterministic evaluations 
based on binary distributions. 

Because each Monte Carlo evaluation can be 
calculated independently, determining stochastic 
robustness is a task well-suited to parallel computation. 
Eigenvalue computation speed is linear in the number of 
processors, and results from separate processors need be 
consolidated only at the final stage of display. 

CONCLUSIONS 

Stochastic robustness offers a rigorous yet 
straightforward alternative to current metrics for control 
system robustness that is simple to compute and is 
unfettered by normally difficult problem statements, such 
as non-Gaussian statistics, products of parameter 
variations, and structured uncertainty. The approach 
answers the question, "How likely is the closed-loop 
system to fail, given limits of parameter uncertainty?'* It 
makes good use of modem computational and graphic 
tools, and it is easily related to practical design 
considerations. The principal difficulty in applying this 
method to controlled systems is that it is computationally 
intensive; however, requirements are well within the 
capabilities of existing computers. The principal 
advantage of the approach is that it is easily implemented, 
and results have direct bearing on engineering objectives. 
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